Attractive Surface Forces due to Liquid

The force between two inert, planar surfaces in a liquid is discussed, focusing on the e ect of the reduction in the uid density between the surfaces. We use a density functional treatment of the uid, the results of which are veri ed by comparison with Grand Canonical Monte Carlo simulations. Our calculations reveal the existence of an attractive force between the surfaces, which at short and intermediate separations is substantially stronger than that expected from the conventional Hamaker or Lifshitz theories. The source of the extra attraction is the interstitial density depression that results from the requirement of chemical equilibrium between the bulk and the con ned uid at all separations. Given inverse sixth power attractions between uid molecules, the force between the inert surfaces asymptotically decays as the inverse third power of the separation. However, at short separations it displays a qualitatively di erent behaviour. We also nd that the attraction increases with temperature, at constant pressure, and decreases with increasing pressure, at constant temperature, contrary to standard Hamaker theory. These results should lead to a more cautious interpretation of some parameters obtained from tting experimental surface force curves to, for example the DLVO expression. Furthermore, in discussions of the molecular mechanism of the hydrophobic attraction, the focus is almost invariably on the short-range directional hydrogen bond, or electrostatic, forces. Our calculations suggest that the density depression mechanism could be as relevant as orientational order. 3 1 INTRODUCTION There are several mechanisms that can generate a long-ranged force between two surfaces (or particles) in solution [1, 2]. For some of these the molecular origin is well understood, and qualitatively correct theories have been developed, which are sometimes also in quantitative agreement with experiments [3, 4, 5, 6]. In other cases there are still unclear points in both theory and experiment. The Surface Force Apparatus (SFA) [5] and the osmotic stress technique [7] as well as the more recent atomic force microscope have been of seminal importance for the development of this area and for the elucidation of di erent contributions to surface forces. The interaction of two charged surfaces, in the SFA achieved by two crossed mica covered cylinders, has been thorougly investigated and it has been veri ed that the DLVO (Dejaguin-Landau-Verwey-Overbeek [3, 4]) theory is qualitatively correct for solutions containing monovalent ions. In particular, the long range part of the repulsive double layer force decays exponentially as predicted by theory. Fitting the experimental results obtained at short separation to the DLVO theory is more troublesome and the determination of, for example, the surface charge density seems less reliable [8, 9]. One reason is of course that the Poisson-Boltzmann (PB) approach inherent in the DLVO theory is a mean eld approximation, which becomes less justi ed at high surface charge density or with multivalent ions where correlations start to play a role. The general trend is that the PB approximation predicts too strong a repulsion leading to an underestimate of the surface charge density. Underestimating the attractive term would of course have a similar e ect on the tted parameters. The interaction between planar hydrophobic surfaces is strongly attractive. However, it has turned out to be very di cult to experimentally determine the range of the force [13, 10, 39, 14, 11, 15, 16, 12]. A survey of the current experimental literature reveals that there is no consensus even on the qualitative features of the distance dependence. This is not the proper place for a detailed discussion of the strengths and weaknesses of 4 di erent experimental studies. Considering that the experimental situation is not settled, this study is made on the tentative premise that the \true" hydrophobic force has a substantially shorter range (<100 A) than some experimental studies have indicated. We will in this study discuss a simple mechanism which gives rise to a strongly attractive force, at short and possibly intermediate separations. Elementary theories of surface forces in liquids assume an incompressible liquid giving a simple relation between surface displacement and deplacement of molecules. This is an obvious simpli cation that can lead to substantial errors. In the present paper we study the consequences for the surface force resulting from depression of the solvent density between the surfaces. This general mechanism is analyzed using a density functional treatment originally formulated by Nordholm and co-workers [20, 21] and which is usually referred to as the "generalized van der Waals theory". The accuracy of this theory is evaluated by comparison with Grand Canonical Monte Carlo simulations from a recent study by Berard et al: [22]. The idea of surface forces being modi ed by density variations induced by the presence of surfaces is not new. Mitchell et al. [23, 24] have previously investigated \structural contributions to Hamaker constants", using integral equation theory. However, their approach is less suited for a study of contributions from a net reduction of the density in the slit, since they used asymptotic forms of the integrals. They also invoked a superposition approximation for the total surface-solvent correlation function, when two surfaces were present. Furthermore, their integral equation theory does not predict a gas-liquid phase transition, which is expected to be of importance as the liquid becomes metastable. Thus, their focus was on e ects of density structure, rather than density reduction. 2 Hydrophobic Interaction The strong tendency for apolar solutes to associate in aqueous solution is formally attributed to the "hydrophobic interaction". For small solutes, the magnitude has been 5 carefully quanti ed both experimentally [25] and via simulations [26, 27, 28, 29], but it has proven di cult to obtain an approximate analytical theory describing the phenomenon. It is clear, however, that the hydrophobic interaction is basically due to the strong cohesion in water, which makes it very unfavourable to generate a cavity where an inert solute can be incorporated. The force between large surfaces, resembling the physical situation in SFA experiments, has also been simulated [30], but such studies are restricted to rather short separations. The surface free energy of an alkane-water interface is approximately 50 mJ=m2 at room temperature. Thus, the total free energy change, i:e:, the integral of the force in bringing two alkane covered surfaces into contact in water should be about 100 mJ=m2. It is, however, not possible to account for such a high surface tension on the basis of van der Waals forces [2]. Consequently, one should not expect that the force between two hydrophobic surfaces is dominated by the van der Waals force, except asymptotically. A simpli ed yet physically appealing picture of what happens in the solvent phase between two inert walls can be obtained by balancing bulk and surface free energy contributions. Consider two identical in nite planar surfaces, S, separated a distance ĥ in a liquid. We denote the interfacial free energy (per unit area) as ̂SL, while the surface free energy relative to the vapour phase is ̂SG. Generally, for solvophobic surfaces, ̂SL > ̂SG, and the surface interactions favour vapour rather than liquid between the walls. Under these conditions there is a possibility of capillary evaporation or cavitation. A simple thermodynamic analysis assuming ideal gas behaviour in the vapour, shows that con ned gas and liquid phases coexist at some separation, ĥcrit, given by ĥcrit = 2( ̂SL ̂SG) P̂B ln h P̂B P̂sati (1) where P̂sat is the saturation pressure and P̂B is the external bulk pressure. For ĥ < ĥcrit, the vapour is the more stable phase. Inserting parameters typical for a hydrocarbon-water surface at P̂B = 1 atm and room temperature, one nds that ĥcrit is approximately 1500 A. This number should be treated with some caution, but it points to a possible source for a long ranged attractive force between two inert surfaces. Prior to evaporation, there 6 will be a lowering of the average density in the slit. Our discussion above suggests that the liquid will be metastable below ĥcrit. Nevertheless, provided the barrier to nucleation is su ciently large, one can treat the metastable uid as being at a stable equilibrium. Metastable conditions are quite common in experimental systems [15] and the free energy functional is able to resolve metastable liquid phases. 3 THEORY 3.1 Model We shall treat a model, in which spherical uid particles are con ned by two surfaces. These surfaces are \hydrophobic" (or more correctly \solvophobic") in that there is no attractive interaction with the uid, only an in nite repulsion at the distance of closest approach. The particle-particle interaction, û(r̂), is given by the Lennard-Jones potential, û(r̂) = 4 [( ̂r )12 ( ̂r )6] (2) The con ned uid is in equilibrium with a bulk at some chemical potential, ̂B. As there is no direct interaction between the surfaces, the force between them is solely due to the action of the intervening uid. The pressure acting on the walls, due to the solvent, can be written as the negative derivative of the free energy per unit area, Ĝs(ĥ), P̂s(ĥ) = dĜs(ĥ) dĥ (3) The solvation pressure will depend on the bulk chemical potential, which in turn is determined by the bulk uid density, n̂B, and temperature, T̂ . The choice of T̂ and n̂B, will also determine the surface tension of the liquid. 7 In the following we will be using scaled variables. In the scaled system, all energies are multiplied by the factor 1= and the temperature is multiplied by kB= , where kB is Boltzmann's constant. Scaled lengths are obtained by dividing by . Physical quantities will be distinguished from their scaled counterparts by the \hat" symbol. 3.2 Generalized van der Waals theory The generalized van der Waals (gdvW) theory [20, 21, 31] has been extensively applied to non-uniform uids. The theory exists in several formulations and the main di erence between them is the length scale over which they are able to resolve uid structure. In the case of a uid con ned by surfaces, the average density pro le will vary most rapidly close to the surface. For a dense hard sphere uid at a hard surface or a dense Lennard-Jones uid close to attractive surfaces, one expects an oscillatory density pro le. In order to resolve oscillatory pro les, one must use a so-called " ne grained" free energy functional [21] in which non-local entropy e ects are included. However, in the case of Lennard-Jones particles close to repulsive walls, one expects signi cantly reduced structure. Further, in this study we are interested in the interaction between surfaces which are separated by distances larger than the typical range of oscillatory structures. Thus it should su ce to use a "coarse grained" free energy functional [20]. These types of functionals are much faster to solve numerically. Using the simplest of the gvdW coarse grained functionals, the free energy per unit area for a particular density pro le n(z), can be written as, Gs[n(z);h] = T Z h 0 dzn(z) ln " 3(1 n(z)) n(z) #+ 1=2 Z h 0 dzn(z) Z h 0 dz0n(z0) (jz z0j) B Z h 0 dzn(z) + PBh (4) The rst term on the right hand side of Eq.(4) is the entropy contribution, as estimated by an approximation of the free volume available to hard spheres with an e ective diameter equal to . The second term is a mean eld approximation to the interaction between uid particles. The nal two terms are the separation dependent bulk contributions. The bulk 8 pressure, PB, and chemical potential, B are related to the bulk density, nB, according to PB = TnB 1 nB 16 n2B 9 (5) B = T ln " 3(1 nB) nB # + T 1 nB 32 nB 9 (6) The laterally integrated pair potential is given by, (z) = 8><>: 4 [ 1 5z10 1 2z4 ] ; z > 1 6 5 ; z 1 (7) The equilibrium density pro le, neq(z), minimizes Gs and thus satis es the following integral equation, n(z) = (1 n(z)) exp " n(z) 1 n(z) 1 R h 0 dz0n(z0) (jz z0j) B T # (8) In this work, Eq.(8) was solved iteratively, until the maximum absolute value of the functional derivative, Gs= n(z), was less than 5 10 6. We used an integration grid of 0.005 in the z-coordinate. We also performed some calculations using half this grid size and no di erence in any property was obtained. A typical calculation was completed in a few minutes on an ordinary work station. For further details we refer to the following work [20, 21, 32] and references therein. The net pressure acting on the surfaces, Ps(h), was evaluated by a discrete numerical di erentiation of Gs with respect to h, i.e by using the scaled quantity version of Eq.(3). The contact theorem for hard walls, [1], Ps(h) = Tneq(0) PB (9) does not apply in the coarse grained versions of the gvdW theory. Interestingly, it does apply for ne grained functionals. Other meaneld approaches have also been used to study non-uniform solvent structure. Of particular relevance this work is the use of the Cahn-Hillard theory to study oil water 9 inerfaces by Turkevich et al. [?, 35]. The Cahn-Hillard theory can be obtained from the gvdW approach by replacing non-local terms with a density gradient approximation. Such an approach will be used later to study asymptotic solutions of the gvdW equations. 3.3 Hamaker Theory In the classical theories of surface forces, due to Lifshitz and Hamaker among others, it is recognized that the interactions between surfaces will be modi ed by the polarizability of an intervening uid [33, 34, 1]. The Hamaker expression for the interaction between the surfaces can be obtained from the gvdW theory by imposing a xed density for the con ned uid, equal to that of the bulk. The free energy functional can then be solved analytically with the following result, Gs;Ham(h) = n2B 1 90h8 1 3h2 (10) The pressure between surfaces in the Hamaker approximation is then, PHam(h) = 2 n2B 2 45h9 1 3h3 (11) The constraint of a constant density pro le for the con ned uid, equal to the bulk value, is a rather drastic approximation, especially close to the surfaces. One of the results of this work is to demonstrate that imposing chemical equilibrium between the con ned liquid and the bulk, leads to a signi cant modi cation of PHam. This was also observed in the simulations of Berard et al., though the emphasis in that work was on the role played by the apparent spinodal separation. Despite the fact that our model only incorporates high frequency responses, it is clear that the e ect we see will impact on the polarizability of the uid at all frequencies, as these all depend on the uid density. 10 4 Results and Discussion We compared systems at di erent temperatures and constant bulk pressure. The bulk density was obtained by solving Eq.(5) at a given temperature. Figure 1 shows the ratio between the gvdW pressure and the Hamaker pressure as a function of separation. As the separation increases, the gvdW and Hamaker results will coincide. This is due to the fact that the average uid density in the centre of the slit will approach the bulk value, with the zones near the walls making an ever diminishing contribution to the free energy of interaction. However, Ps is much more attractive than PHam at short and intermediate separations. The two pressures can di er by more than an order of magnitude at short separations, and the discrepency can still be as large as 50 per cent at a separation of 40 molecular diameters. One can interpret this additional attractive force (beyond the Hamaker component) as the contribution due to uid structure. That is, it arises due to the deviation of the con ned uid density away from the bulk value. The Hamaker pressure, PHam, decays as the inverse third power of separation, but it appears that the structural component of the force decays more quickly, at least at short separations. By studying Ps PHam we found that the extra contribution decays algebraically, with the leading long-range term varying as h 4. We present below a perturbation expansion which con rms these ndings. At distances less than the minimum separations depicted on the gures, no stable liquid phase solution to the gvdW equations could be found. At these minimum separations the liquid becomes spinodal. Recently Berard et al., suggested that the spinodal separation could in uence the liquid correlation length (even at much larger separations) and that this may account for the long-ranged hydrophobic interaction. We only note here that any such e ect (regardless of its importance) will be implicitly included in the gvdW theory. The validity of Eq.(1) was roughly estimated by inserting gvdW values of the required quantities and comparing the predicted critical separation, as found from this expression and gvdW calculations, respectively. Not surprisingly, the agreement is best at high 11 0 20 40 h 0 2 4 6 8 10 P S / P H am Figure 1: The ratio between the gvdW pressure, Ps from Eq. (9) and the Hamaker pressure, Eq. (11), as a function of separation . Temperature and distance are given in reduced units and PB = PB;ref = 0:11074 0:00003. Three di erent temperatures are shown: T = 0:9 as a dot-dashed, T = 1:1 as a dashed and T = 1:3 as a solid line. temperatures, where the bulk pressure is closer to the saturation pressure. At T = 0.9, Eq.(1) and gvdW calculations then predicts hcrit = 8 and 21, respectively. At T = 1.3, the agreement is much better, with hcrit = 19 and 27, respectively. A special feature of the observed attraction is that, for the parameter range investigated, it appears to get larger with increasing temperature. This is opposite to the behaviour of the Hamaker component. Increasing the temperature at constant bulk pressure leads to a reduction of the bulk density. Eq.(11) indicates that this leads to a decrease in PHam. The additional structural component gives a much more attractive force at higher temperatures, at short to intermediate separations. At larger separations, the structural component has decayed signi cantly and the total force approaches the Hamaker limit, see Figure 2. In SFA experiments, which use crossed cylindrical surfaces, the measured force is the net surface free energy of interaction, Gs(h) 2 SL, between at surfaces (assuming the Derjaguin approximation). In the numerical solution of the gvdW equations, we chose to set Gs(h) 2 SL equal to zero at h = 200, which is large enough to give a very good estimate to the limiting surface tension. The resulting free energy curves are given in Figure 3. Again 12 0 20 40 60 80 h −6 −4 −2 0 P *h 3 Figure 2: The gvdW pressure and the Hamaker pressure multiplied by h3 as a function of h at PB = PB;ref , for three di erent reduced temperatures: T = 0:9 as a dot-dashed, T = 1:1 as a dashed and T = 1:3 as a solid line. The same symbols have been used for the thin horizontal lines representing h3PHam. 0 20 40 h −0.02 −0.01 0.00 G S ( h) − 2 γ S L Figure 3: Free energies of interaction, Gs(h) 2 SL, as obtained using the gvdW theory (thick lines) and Hamaker approach (thin lines), respectively. The results for three different reduced temperatures are given: T = 0:9 as a dot-dashed, T = 1:1 as a dashed and T = 1:3 as a solid line. we see that, as the temperature is increased at constant bulk pressure, the gvdW theory predicts a stronger attraction at shorter separation, in contrast to the Hamaker result. A 13 study of the density pro les obtained from the gvdW theory, gives an indication of the mechanism behind this behaviour. Figure 4 displays density pro les, obtained at a xed separation. It clearly shows a density depletion close to the walls and how, at higher temperatures, the e ect of the walls propagates further out. At rst sight this result seems counterintuitive. A reduction in the importance of the interactions between particles, as would occur upon raising the temperature, should reduce the depletion e ect. Indeed, if the bulk density were to remain xed and the temperature were to become in nite, there should be no depletion. As we shall see below, however, it is important to keep in mind that, in the systems studied above, the bulk pressure is xed not the bulk density. This boundary condition is more consistent with experiments. With the pressure xed, increasing the temperature increases the correlation length in the liquid, thus leading to an extension of the density depleted region towards the middle of the pore. Figure 5 shows how the average density varies with slit width and temperature. Again we see that the e ect of the wall increases with temperature. It is worthwhile pointing out that we also repeated some density pro le calculations using ne grained free energy functionals. We found very little di erence in the results, provided the comparison was made at the same T=Tc, where Tc is the critical temperature as predicted by the particular version of the theory. This is not surprising, given that the uid con ned by repulsive surfaces does not possess oscillatory structure. In order to better understand the temperature e ect on the uid density pro le, and consequently the free energy of interaction, it is useful to consider a perturbation approximation to the gvdW theory. In particular, we expand the functional Gs[n(z);h] in Eq.(4) about the bulk density, to second order in n(z). This should be a reasonable approximation for the uid density pro le far from the walls for large h. We show in the Appendix that the leading order expression for the density pro le around the centre of the slit is 14 0 4 8 12 z 0.0 0.2 0.4 0.6 0.8 1.0 n( z) Figure 4: Density pro les for three di erent temperatures at a reduced separation of H = 12 and PB = PB;ref : T = 0:9 as a dot-dashed, T = 1:1 as a dashed and T = 1:3 as a solid line. The corresponding bulk densities are given as thin horizontal lines. 0 20 40 60 80 100 h 0.7 0.8 0.9 1.0 n( h) / n B Figure 5: The averaged reduced density as a function of separation for three di erent temperatures, at PB = PB;ref : T = 0:9 as a dot-dashed, T = 1:1 as a dashe d and T = 1:3 as a solid line. given by n(z) nB 1 (z) (12) 15 where (z) is given by Eq.(20) and, = @ B @nB!T = T nB + T 1 nB + T (1 nB)2 32 9 (13) The function (z) accounts for the e ect of the walls on the density pro le. It is positive and decreases monotonically away from the walls, depending on the inverse third power of the distance from the walls. The magnitude of its in uence is determined by 1= . The correlation length of the bulk uid is proportional 1=2. Thus altering conditions so as to increase the correlation length, will increase the e ect of the walls. For a liquid below the critical temperature, increasing the temperature at constant pressure will move it closer to the bulk liquid/vapour phase boundary. The gvdW equation of state predicts that this leads to an increase in the correlation length, as the liquid phase becomes in some sense \less stable". This explains why, in the example calculations above, the e ect of the walls propagates further out at higher temperature. On the other hand, if the bulk uid is a gas, increasing the temperature at constant pressure will reduce the correlation length, as the uid moves away from the bulk phase boundary. For the same reason, if the temperature is increased at xed density in either a bulk gas or liquid, the correlation length decreases. Note that at critical or spinodal conditions in the bulk uid, one has = 0, and the correlation length becomes in nite. There have been some attempts to experimentally measure the width of uid interfaces and it may be interesting to compare with these. We have attempted to calculate the width of the interface layer against a single surface using the formula [35], = @ ̂ @P̂B (14) For the pentane/air interface the experimental value [36] is = 2.2 A. From the gvdW theory we obtained a value of 1.6 A. The calculations employed Lennard-Jones parameters derived by Rodriguez and Freire [37]. From the perturbation expansion of the gvdW functional the free energy is given by, Gs(h) = Gs(nB; h) + 1 2 Z h 0 dz neq(z) (z) (15) 16 where Gs(nB; h) is the free energy obtained by assuming a constant uid density, equal to the bulk value. The separation dependent part of Gs(nB; h) is the Hamaker interaction, Eq.(10). The second term in Eq.(15) corrects for uid structure. We can approximate it with the leading order term for the density, to give Gs(h) = Gs(nB; h) 1 2 Z h 0 dz (z)2 (16) The structural correction term contains surface and separation dependent components. The surface component lowers the surface tension term obtained by assuming a constant liquid density pro le. One can show that the separation dependent part, asymptotically varies as h 3 and thus decays faster than the Hamaker component, which goes as h 2. As mentioned earlier, such a decay was indeed obtained from our gvdW calculations, but at separations less than about 40-50 molecular diameters, higher order terms contribute signi cantly. The structural part of the free energy of interaction also increases with the correlation length. This increase in attraction can more than compensate the reduction in the Hamaker component, as was observed in the examples above. Note that Mitchell et al. [23], using an entirely di erent approach, obtained qualitatively identical results, at large separations. According to the above discussion, lowering the bulk pressure at constant temperature will also move the thermodynamic state of the liquid closer to the phase boundary and may thus increase the attraction at short to moderate separations. This is con rmed in Figure 6, where a fty per cent reduction of PB leads to a doubling of the attraction. An important test of the generalized van der Waals theory is a direct comparison with Monte Carlo or Molecular Dynamics simulations. Fortunately, a Grand Canonical Monte Carlo study of a Lennard-Jones liquid between two hard walls has recently been published by Berard et al. [22]. We performed test calculations using the gvdW theory with precisely the same parameters and were able to reproduce the MC results within the simulated error bars, Figure 8. The excellent agreement (also found for the density pro les) suggests that one can use the gvdW theory with a great deal of con dence in these types of systems. It 17 0 20 40 60 h 1 3 5 P S / P H am Figure 6: The ratio between the gvdW pressure, Ps from Eq. (9) and the Hamaker pressure, Eq. (11), as a function of separation , at T = 1:1. Three di erent bulk pressures are shown:PB = 2PB;ref as a dot-dashed, PB = PB;ref as a dashed and PB = PB;ref=2 as a solid line 0 20 40 60 80 100 h −3 −1 P * h 3 Figure 7: The gvdW pressure and the Hamaker pressure multiplied by h3 as a function of h for three di erent reduced pressures: PB = 2PB;ref as a dot-dashed, PB = PB;ref as a dashed and PB = PB;ref=2 as a solid line. The same symbols have been used for the thin horizontal lines representing h3PHam. should be pointed out that the computation time for the simulations was several orders of magnitude greater than required to solve the gvdW equations. 18 4 5 6 7 8 9 h −0.2 −0.1 0.0 P S Figure 8: GvdW (solid line) and GCMC (diamonds) [22] data for a Lenna rd-Jones liquid con ned between two hard walls, at T = 1:0 and PB = 0:44. In the GCMC simul ations, the bulk pressure (and thus the net pressure) is only determined within 0:05. It would be of considerable interest to attempt to model liquid water, and to estimate the e ect that the density depression mechanism would have in that system. However, the gvdW theory, at least in any of its simpler versions, is essentially limited to studies of spherical molecules. The approach we take here is to use an e ective Lennard-Jones uid, with some properties similar to that of water. We considered two possible e ective models. In the rst, we used gvdW results at the same T=Tc and P=Psat as bulk water at room temperature and atmospheric pressure, Tc and Psat being the critical temperature and saturation pressure, respectively. In the second model, the interfacial surface tension, , and the bulk isothermal compressibility, T , as given by the gvdW theory were chosen to be similar to those of water. The results using these two models are given in Figure 9. We stress that one could quite validly propose the use of some other approach e ective model, given the non-uniqueness of such a mapping to the spherical uid. However, the di erences that can be expected are, to some extent, indicated in Figure 9. We also calculated the value of for these models, using Eq. (14) and obtained the values 2.5 and 2.9 A, repectively (assuming a molecular diameter of 3.1 A). The experimental value for water at the liquid-gas interface interface is 7.0 A [36]. A general feature of both models is 19